MacFormat 1994 August/August CD.bin/Shareware/Education/numericalmethods Folder/chap_10/crnich.m

Text File | 1994-06-05 | 1018 b | 43 lines | [MATF/MATL]

function U = crnich(f,g1,g2,a,b,c,n,m)
% U = crnich(f,g1,g2,a,b,c,n,m)
% Crank-Nicholson solution to the heat equation.
% This subroutine uses trisys.m to solve a tri-diagonal system.
% f  is a function, input.
% g1 is a function, input.
% g2 is a function, input.
% a  is the width of [0 a], input.
% b  is the width of [0 b], input.
% c  is the constant in the heat equation, input.
% n  is the number of grid points over [0 a], input.
% m  is the number of grid points over [0 b], input.
% U  is the solution matrix, output.
h = a/(n-1);
k = b/(m-1);
r = c^2*k/h^2;
s1 = 2 + 2/r;
s2 = 2/r - 2;
U = zeros(n,m);
for j=1:m,
     U(1,j) = feval(g1,k*(j-1));
     U(n,j) = feval(g2,k*(j-1));
end
for i=2:(n-1),
     U(i,1) = feval(f,h*(i-1));
end
Vd = s1*ones(1,n);
Vd(1) = 1;
Vd(n) = 1;
Va = -ones(1,n-1);
Va(n-1) = 0;
Vc = -ones(1,n-1);
Vc(1) = 0;
Vb(1) = feval(g1,k*0);
Vb(n) = feval(g2,k*0);
for j=2:m,
  for i=2:(n-1),
       Vb(i) = U(i-1,j-1) + U(i+1,j-1) + s2*U(i,j-1);
  end
     X = trisys(Va,Vd,Vc,Vb);
     U(1:n,j) = X';
end